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ABSTRACT, 

A numerical model which utilizes the isobaric vorticity equation 
is developed and applied to cumulus -s ca le data. The model, together 
with a modified version of the cumulus convection model of Weinstein 
and Davis, is applied to data obtained from the National Severe Storms 
Laboratory in Norman, Oklahoma. The calculations yield real time pre- 
dictions for height, temperature and relative humidity at seven pres- 
sure levels, which are then used as input to the cumulus convection 
model to obtain vertical profiles of various parameters at specified 
grid points . 

Some results of the calculations are presented along with sug- 
gestions for further testing and improvement. The results indicate 
that further modifications to the approach used are necessary in order 
to provide more accurate forecasts. Values of the individual terms 
in the vorticity equation are presented as computed from the observed 



mesoscale data. 
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1. INTRODUCTION 



The National Severe Storms Laboratory (NSSL) located at Norman, 
Oklahoma, has been gathering vast amounts of data during the thunder- 
storm season for the past several years. Many studies have been 
conducted with the aid of this data (for example see Barnes, 1968, 
Newton and Fankhauser, 1964 and Hammond, 1967), however, after an 
extensive literature search, it appears that no research has been 
done, or at least reported, involving a numerical modelling approach 
whereby the distributions of various parameters were analyzed and their 
values numerically forecast for a future time. Such an approach, if 
successful, could give distinct advantages in predicting, for example, 
where hail might be expected, or where the most favorable conditions 
for tornadoes might be located. In addition, during severe convective 
activity radiosonde data often terminates near mid- tropos pheric levels. 
Consequently, a model which could reasonably predict the sounding 
above termination level, when given a complete sounding from the 
previous observation time, would aid greatly in later research efforts. 

Figure 1 shows the region used for this study with a height-field 
analysis and the reported winds at the nine radiosonde stations in 
the NSSL upper-air network (see Appendix A for a description of the 
NSSL upper-air network). Assuming the data to be free of errors, 
the non-geos trophic nature of the flow on this scale is immediately 
evident, and the figure emphasizes the need for intensive research 
into the dynamic mechanisms involved under such circumstances. In 
an effort to evaluate the possibilities of numerical weather pre- 
diction for such small-scale data, this paper uses a numerical 
prediction model based on the complete isobaric vorticity equation. 
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Figure 1: Height analysis with wind directions and speeds 
at nine radiosonde stations of NSSL network for the 500-mb 
1700Z on 30 May 1967 which illustrates the non-geos t roph i c 
the da ta . 



(m sec ; 
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The resulting forecasts are then used as input to a cumulus convection 
model developed by Weinstein and Davis‘S (1967) to provide an aid to 
severe storm forecasting. 

The parameterized cumulus convection model by Weinstein and Davis 
(1967) calculates cloud top, profiles of vertical velocity, tempera- 
ture excess (that is, the temperature of a parcel of air in the cloud 
minus the temperature of a parcel of air in the environment), hydro- 
meteor water content and updraft radius of the cloud as well as rain- 
fall amount and duration. The principles of the above calculations 
are based on the entrainment principles of Stommel (1947) and the 
results of Kessler (1967). With a standard radiosonde sounding at 
mandatory and significant levels beginning at cloud base, the W.D. 
model transforms the sounding to one which is spaced in equal incre- 
ments of height. After the user specifies certain boundary conditions 
such as ice nucleation temperature, radius of the base of the cloud 
and conversion and collection rates (Kessler, 1967), the model 
specifies cloud extent. One use of this model is to compare the 
effect of seeding a given cloud or allowing it to develop naturally. 
This is accomplished by varying the ice nucleation temperature. 

The numerical model presented in this paper predicts heights, 
temperatures and relative humidities at individual grid points for 
seven isobaric levels (950 mb, 850 mb, 700 mb, 500 mb, 300 mb, 200 
mb, and 150 mb) so that predicted information can be used as input to 
a modified subroutine form of the W.D. model to yield the vertical 

^The model developed by Weinstein and Davis has been referred to as 
the Pennsylvania State model. In this paper, it will be abbreviated 
as the W.D. model. 
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profiles previously discussed, under the assumption that cloud base 
is at 850 mb. However, the present version, as described in this 
paper, does not provide for vertical coupling of the seven isobaric 
levels after the first time step. 

Predictions were made at 1.5-hour intervals utilizing the IBM 
OS/360 (MVT versions 15, 16) computer. The 1.5-hour prediction 
interval was selected to correspond to the times when verification 
data were available from the NSSL radiosonde network. 

The grid used for the computations is a 24 x 24 X-Y grid oriented 
000-180, 090-270, having a mean latitude of 35N. The grid distance 
is 5 nautical miles with the grid arranged so that the nine upper- 
air stations are either located exactly on a grid point or are so 
close to one that they may be approximated to be at the nearest grid 
point. The farthest distance away from a grid point for a particular 
station was 2 nm. They were assumed to be located at the nearest 
grid point to eliminate the need for an interpolation scheme in the 
objective analysis used for initialization of the data as described 
in Section 4. 
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2 . MODIFICATIONS TO THE W.D. MODE L 



While working with and studying the W.D. model, it became apparent 
that several modifications could be made to it that would provide more 
accuracy. Although most of the modifications were of a minor nature, 
three of the modifications are considered to be significant improve- 
ments over the original model. These three modifications will be 
presented in detail while the remainder of the changes will be sum- 
marized at the end of this section since they are considered to be 
se If -explanatory . 

Kessler (1967) proposed the equation 



to model the conversion of cloud water to hydrometeor water. Here, 



and is zero for m<a . The quantity "a" is defined as a threshold 
value and refers to the cloud water content. If the cloud water 
content is greater than this threshold value, then conversion of 
cloud water to hydrometeor water will occur. If cloud water amounts 
less than the threshold value are present, then conversion of cloud 
water to hydrometeor water will not take place. This equation is 
incorporated into the W.D. model (Weinstein and Davis, 1967) in the 



(1967) indicated that measurements have been performed which indicate 

-3 

that cloud water amounts greater than 1 gm m are usually associated 




(2 -la) 



k^ is an empirical constant which is equal to .001 sec when m>a 



form 




(2-lb) 



where k^ is defined in a similar manner as it was in (2-la). Kessler 
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with the production of precipitation. Furthermore, he arbitrarily 

-3 

suggested the value of 0.5 gm m for the threshold value. Weinstein 
and Davis (1967) chose a* = 0.5 gm Kg ^ for their threshold value. 

It is proposed to consider this threshold value to be a function 
of density. This approach thereby associates a different threshold 
value with each level by making it a function of height. Using this 
approach, the conversion of cloud water content to precipitation in 
the W.D. model has been modified to 

dQ h 

— - k 1 «c- /p > <2 - 2) 

where p = air density. By defining as it was in (2-la), conversion 
of cloud water to precipitation is modelled to occur when (Q^ - — ) > 0. 

Computations have been made to determine the variation of the 
quantity a/p. These values show a maximum of approximately .00137 
gm gm ^ in the upper levels and a minimum value of .00076 gm gm ^ 
near cloud base. 

The second significant modification is the inclusion of a linear 
variation of the latent heats of vaporization (L^) and fusion (L^) 
with temperature. Since 

L - L + L (2-3) 

s v r 

and L g may be considered constant while is a function of temperature, 
then L is also a function of temperature. The values of (List, 

1951) are 48.6 cal gm ^ at -50C and 79.7 cal gm ^ at 0C . The latent 
heat of vaporization varies from 629.3 cal gm ^ at -50C to 597.3 
cal gm ^ at 0C . Since L g varies by only 1 cal gm ^ from 0C to -40C , 
its value was taken as constant and equal to 678 cal gm The 
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variation of only one of the latent heats on the right side of equation 
(2-3) needs to be computed since the value of the second one can be 
obtained by subtraction. 

In order to approximate the variation of the latent heat of vapor- 
ization as a linear function of temperature, it was necessary to deter- 
mine an optimum value for the difference between the specific heats of 
water and ice. Integrating the formula 



dL 

- -r = C - c 

dT pv w 

and applying the mean value theorem for integrals yields 



L = L + (C - C ) (T-T ) 
V V pv w' v o 

o r 



(2-4) 



(2-5) 



Utilizing the values of (List, 1951) for various temperatures (T) 



and with T = -40C, the value of (C - C ) that satisfies equation 
o pv w 

(2-5) was computed for each temperature at 10 degree intervals from 



-40C to +40C. A plot of (C - C ) versus temperature is shown in 

pv w r 

Figure 2. Based upon the results it was decided to select the value 



of -.62340 as an "optimum" value for (C - C ) and consider it to be 

r pv w 

a constant. This value is considered to be a representative value 
between -40C and 0C , and corresponds exactly to the value at -10C. 

The final significant modification involves the values used for 
the constants when computing saturation vapor pressure (Weinstein 
and Davis, 1967). Integrating the C laus ius -C lapeyron equation one 
obtains the following for the saturation vapor pressure over water: 



e [L - T (C - C )] € (C - C ) In T 

_ VO o pv w + pv w 

n e s RT R 



e [L - T (C - C )] e (c - C )ln T 



VO 



o pv w 

RT 






w 



R 



+ In e so . (2-6) 
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Simplifying yie Ids 



In e = 
s 



- — - B' In T + C' 



where 



A' = 



c [L -T (C - C )] 



VO o pv 



w 



= 6958. 9262 



e (C - C ) 

, _ pv w 



B = 



R 



= 5.65567 



and 



C' - § |> - (C - C ) (1+ln T )~| + In e = 59.01383 
R LT pv w o J so 



when 

and 



lo - 1 



T = 233. 16K, (C - C ) = -0.62340 cal gm K 
o pv w 

L =621.7 cal gm ^ . 
vo 



Similarly, when freezing occurs one obtains 

-A M 



lne s=t + C " 



(2-7) 



where 



eL 



A" = 



R 



= 6151.0205 



A" -1 

and C n = — + In e = 24.3277 when T = 233 . 16K and L = 678 cal gm 
T so o s 

o 

The above equations yield the saturation vapor pressure in millibars. 
Obviously the number of significant digits shown above exceeds the 
accuracy of some of the assumptions made to arrive at the values shown. 
Rounding to five significant figures could be done with no loss of 
accuracy, however they have been retained here as they were computed 
on the IBM OS/360 computer. 

There were several minor changes made to the W.D. model. These 

included changing the value being used for the acceleration of gravity 
-2 

from 9.87 m sec to a more exact form obtained by computing the value 
from the formula (List, 1951) expressing g as a function of latitude, 
altering the values used for the constants for conversion from degrees 
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Celsius to degrees Kelvin, conversion from meters to feet, and the 
value of the gas constant for dry air. In addition, the section of 
the program dealing with corrections due to the wind shear were com- 
pletely eliminated which consequently allowed a subroutine (subroutine 
CHEQ) to be eliminated. Other minor changes included removing factors 
of 10 so that all pressure values are carried in units of millibars, 
introducing symbols wherever possible to simplify future modification 
and testing, changing the estimate of the updraft velocity at cloud 
base from 2.0 to 0.5 m sec \ and modifying certain format statements 
to improve the headings on the output. Finally, many soundings were 
run through the final modified version of the W.D. model to ensure that 
all changes had been properly programmed and to ensure that the results 
from the model were reasonable. 
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3. MATHEMATICAL DEVELOPMENT OF A FORECAST MODEL 



This section describes the model used for this research. Although 
it is not consistent in many respects, it did provide a means by which 
forecasts of height at 1.5-hour intervals could be obtained for input 
to the W.D. model as well as provide a method by which the values of 
the various terms in the prediction equation could be ascertained. A 
description at the end of this section describes an attempt to set up 
the prediction equation in a consistent manner after the original 
model was running properly, however, the more consistent method failed 
to provide reasonable forecasts. 

A. THE VORTICITY EQUATION 

Following the procedures given by Haltiner and Martin (1957), the 
isobaric vorticity equation may be written as 

&«♦*>-- «♦**■«-♦(§£-&£>. (3- la) 

Since the entire grid used for this research covers only 2 degrees 
of latitude, the Coriolis parameter can be considered constant, 
taking on the value it would have at the mean latitude of the grid. 

Using this approximation and expanding the left side of (3-la) yields 

I=.( u | +W . vc+ (c+ f m)v . W + (3- lb) 

In order to use the W.D. model at 1.5-hour intervals, height of the 
base of the clouds is needed from which computations can be started. 

To obtain these heights, (3-lb) will be used as a height forecasting 
equation analogous to methods applied to large scale motions in the 
atmosphere by applying the balance equation. An alternative method 
without using the balance equation for initialization was tried and 
is discussed at the end of this section. 
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Since the application of (3-lb) to height forecasting will be done 

to obtain information regarding the orders of magnitude of the various 

terms and to obtain heights which can be used as cloud base information 

for input to the W.D. model, no attempt was made to set up a consistent 

set of equations. To obtain an equation which could be used for height 

forecasting and evaluation of the various terms with mesoscale data, 

the relative vorticity on the left side of (3-lb) and in the diver- 

2 

gence term was replaced by V i|f. In addition, only the rotational wind 
components were used for the horizontal advection term, that is, 



di|f 

dy 



u = - and v = 



(3-2) 



This was done because the present version has no provision for pre- 
dicting the divergent component of the wind ( \V ). Future modifications 

X 

should provide the means by which \V can be used in addition to fore- 
casts of the rotational component (\V^). The vorticity in the vertical 
advection term as well as divergence and the u and v components of the 
wind in the twisting term were computed based on the observed winds 
at the initial time and held constant for 1.5 hours. Keeping values 
constant in the terms greatly simplifies programming considerations 
since 3-minute time steps were used and to update these values at each 
time step would involve working with each level for only one time step 
before proceding to the next level. Unfortunately, holding Q in the 
vertical advection term and the u and v components of the wind in the 
twisting term constant for 1.5 hours causes the seven levels to lack 
vertical coupling. Future modifications should attempt to correct this 
deficiency. After making these substitutions and assumptions, the 
resulting prediction equation becomes 

SF< v2 *> - - C4 - m + 1 1 + <" 2 * + 
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Assuming continuity in if , the term on the left side of (3-3) may 
be written as *^, and the equation then becomes a Poisson-type 
equation which can be solved by relaxation procedures as described by 
Haltiner and Martin (1957). The forecasting procedure, being one of 
successive iterations, will thus allow all terms involving in (3-3) 
to vary at each time step. 

To obtain a stream field for (3-3), consider the balance equation 
(Charney, 1955) : 






(3-4) 



Since the Coriolis parameter is to be considered constant over the 
grid, the second term of (3-4) drops out. 

As discussed by Charney (1955) and Bolin (1955), equation (3-4) 



will maintain its elliptic character only if 



V $ f 

— + 2 >0 ’ 



(3-5) 



that is, only if the geostrophic relative vorticity is greater than 
-f/2. As shown in Appendix B, heights would have to be read at least 
to the nearest centimeter on this size grid in order to satisfy the 
requirements of (3-5) and still allow negative geostrophic relative 
vorticities. Consequently, (3-4) was forced to be elliptic by 
writing the non-linear terms in the form 

° ^) m _ dv 

yj Sx Sy 5 



and 



.SxSy> 



cMjf _ Sv Su 

5x 2 8y 2 = 



(3-6a) 

(3-6b) 



Making this substitution into (3-4) and taking into account that 
the Coriolis parameter is to be considered a constant, (3-4) becomes 

V 2 t = §- V 2 z - J- J (u, v) , (3-7) 



m 



m 
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where J represents the Jacobian operator, and u and v are the components 
of the observed winds. Equation (3-7) is then an elliptic equation as 
it stands and does not require that the geostrophic relative vorticity 
be greater than -f/2. The t field can be initialized by the approxi- 
mation and then equation (3-7) can be relaxed to convergence 

m 

before using (3-3) to predict the fields at a later time. Similarly, 
once the | fields have been forecast ahead, the height fields can be 

f m* 

recovered by making the approximation z = ^ and then utilizing (3-7) 

in the relaxation procedure until convergence is attained. 

Convergence was defined as occurring when (3-7) balanced within 

0.5 meters. When converting from height fields to stream fields 

2 -1 

this requires an epsilon of approximately 58612 m sec , whereas 

epsilon becomes 0.5 m when converting from stream fields to height 

fields. These values arise from the approximation i|/ = (g/f )z. In 

m 

retrospect, the above convergence criteria may be too large for this 
scale of motion. In the future, more stringent criteria (such as 0.1 
m) should be invoked. 

As mentioned at the beginning of the section, attempts were made to 

correct the deficiencies of the prediction equation. This involved 

2 

utilizing the exact relation V i|f = Q as the method of obtaining the 

initial field instead of using the balance equation and to substitute 
2 

V for Q in the horizontal advection term. In addition, the conver- 

2 -1 

gence criteria was changed to 10000 m sec , in units of \|f , when con- 
verting from Q to . 

The balance equation was used to recover height fields from the 

forecast fields of i(f with convergence criteria of 0.1 m. Two methods 

of initializing the i|f field were tried; they were ^ and ^ 9 

m m 
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where z represents the standard atmosphere value of the height of the 
respective isobaric surfaces. When using z for the initialization, 
computations were carried out without difficulty, however the result- 
ing forecast height fields were very flat as might have been anticipated. 
However, the forecast fields were in error by as much as 300 m at some 
levels. When using z for initialization, instability resulted in the 
horizontal advection term which caused an exponential growth of this 

term to values too large for the computer to handle after only 12-15 

2 

time steps (36-45 minutes ahead in time). The V part of the hori- 
zontal advection term is believed to be the cause of the instability 
but this was not investigated in detail. However, one attempt to 
investigate the instability was made by shortening the time step to 
one minute. This delayed the instability somewhat which would indicate 
the problems were caused by computational instability. 

B. VERTICAL MOTION COMPUTATIONS 

In order to solve (3-3), vertical motion values are required for 
the vertical advection of relative vorticity and for the twisting term. 
Vertical motions for each of the seven levels were computed by means 
of the continuity equation in the form 



where the bar symbol denotes a vertical average through the layer. 
This method of computation utilizes the measured wind fields and does 
not allow the o> values to be updated for forecast times since the 
present model has no provisions for obtaining the divergent component 
of the wind. The vertical boundary conditions imposed were 




(3-8) 



05 



1000 



= 0.0 



(3-9a) 



and 



150 



0.0 



(3-9b) 
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150 

which also imply that V* W == - Z V* \V. 

1000 p = 950 



Vertical motion computations were made by working downward from 
150 mb through 300 mb and upward from 1000 mb through 700 mb. Vertical 
motion at 500 mb was then computed as the average of the values at 300 
mb and 700 mb. It is recognized that using a backward or forward finite 
difference form for could at times lead to computational difficulties, 
however the method used here does allow the vertical motion at a given 
level to be affected by not only the divergence at that level and the 
level above, but also takes into consideration the divergence field at 
the level below. 

Computations for the various levels were as follows (symbols such 
as W150 and W950 were chosen to conform to the symbols used in the 
computer program (subroutine VERTMO) shown in Appendix C) : 



(1) 150 mb 
W150=0 . 0 

(2) 200 mb (Ap>0) 
Wl75=Wl50-ApVMV 200 _ 15() 
W225=Wl50-ApVMV 300 . 150 
W200=(W17 5+W225)/2 

(3) 300 mb (Ap>0) 
W250=W200 -Ap^W 300 . 200 
W350=W200-A P ^W 500 _ 200 
W3 00= (W2 50+W3 50 ) /2 

(4) 1000 mb 



W 1000=0.0 
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(5) 950 mb (Ap<0) 

W925=W1000-A P ^W 850 . 1000 

W975=W1000-Ap^ 950 _ 1000 

W950=(W925+W975)/2 

(6) 850 mb (Ap<0) 

W825=W950-Ap^ 700 _ 950 

W900=W950-A P ^1v 850 _ 950 

W850=(2(W825)+W900)/3 

(7) 700 mb (Ap<0) 

W675=W850-Ap^lV5 00 „ 850 
W725=W850.Ap^ 500 „ 950 
W700=(W67 5+W725) /2 

(8) 500 mb 

W500= (W3 00+W7 00 ) /2 

The values obtained for vertical motion by this method are discussed 
in Section 5. 

C. TEMPERATURE AND RELATIVE HUMIDITY FORECAST SCHEME 

In order to utilize the W.D. model as an aid in severe storm fore- 
casting, forecasted values of temperature and relative humidity are 
necessary in sounding form as input to the model. 

Although a model designed for operational use should definitely 
be internally consistent, such as the vertical consistency of tempera- 
ture fields with the geopotential fields, this study involved deter- 
mining some of the problems that would develop when applying numerical 
forecasting techniques to mesoscale data. Consequently, the scheme 
for obtaining temperature forecasts which follows is not dynamically 
consistent with what has preceded, but it does provide a method of 



29 



obtaining a tempera ture- fie Id forecast which could be used as input 
to the W.D. model. As discussed in Section 6, future studies should 
improve on this method. 

Four different methods of forecasting temperature were tested and 
the results were then compared to verifying data. The method which 
verified best was selected as the forecast method for temperature, 
which was simply 



where the values for the wind components were taken as the measured 
winds at observation time and were not altered during the forecast 
interval. 

Similar testing of four prediction methods was also done for 
relative humidity. The results led to 



as the prediction equation for relative humidity where the wind 
components were again taken as the measured winds at observation 
time and were held constant for the forecast interval. 

Temperature and relative humidity prediction calculations by 
equations (3-10) and (3-11) were performed in subroutine PROG, a 
listing of which may be found in Appendix C. Results of the com- 
putations, when compared to verifying data, are found in Section 5. 

D. BOUNDARY CONDITIONS 

To initialize the i|r field, the approximation \|r ~ gz/f^ was used 
as discussed in Section 3A. When the relaxation was performed, the 
outside row or column of the grid could not be modified due to finite 




(3-10) 




= - l/2\V-VR h 



(3-11) 
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differencing. As a result, the boundary values of were held at the 
value obtained by the initialization. The prognostic equation (3-3) 
was utilized from x = 3 , 22 and y = 3, 22. After initializing the 
^ field to zero, it was then relaxed for the first time step. This 
relaxation was performed from x =* 3 , 22 and y = 3, 22 as mentioned 
above and the outer two rows or columns were maintained at their 
original values of zero. By this means, actual forecast values of 
| were obtained only from x = 3 , 22 and y = 3, 22. After the first 
time step, the previous relaxed values for were used as the initial 
estimate for the new time step. 

After obtaining the 1.5-hour forecast for the ^ field, the con- 
version back to z fields by means of (3-7) was accomplished by making 
the initial approximation z ~ f^ ty/g only for the points where the ^ 
field had been forecast. This then leaves the outer two rows or 
columns without forecasted height values. Before relaxing the z 
field, however, the values at these boundary points were set by 
performing a linear extrapolation from the fourth row or column 
outward as follows: 



4 . 




3 . 




2 . 


z 2 = z 3 - (z 4 - z 3 ) 


r— 1 


Z 1 = z 2 ' (z 3 ' z 2 } 



Once this was accomplished at all of the grid points in the two 
boundary rows or columns, the z field was relaxed by means of the 
balance equation from x ~ 2, 23 and y = 2, 23. 

Horizontal boundary values were not a problem with the vertical 
motion fields since they were computed from x = 2, 23 and y = 2, 23. 
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As mentioned above, the prognostic equation was utilized from x = 3, 22 



and y = 3, 22 so that there were computed values of omega at x = 2, 
x = 23, y = 2 and y = 23 to use for the terms and As stated 

in (3-9), vertical boundary conditions for a) were that the vertical 
motion would vanish at 150 mb and 1000 mb, which then implies that 
the divergence at 1000 mb is equal to the negative of the sum of the 
divergences from 950 mb to 150 mb. 

As will be discussed in Section 4, the initial data fields of 
height, temperature, relative humidity, wind direction and wind 
speed were analyzed on the computer by an objective scheme. The 
analysis method included finite differencing methods such that the 
outer row or column was not analyzed. These boundary values were 
set to the same value as the next interior row or column. 

The same linear extrapolation technique that was used for the 
height fields was also used for the boundary conditions of the fore- 
casted values of temperature and relative humidity, but it was only 
necessary to extrapolate from the second row or column out to the 
first row or column. 



E. TIME STEPS 

In the forecasting of the new fields by (3-3), a 3-minute time 
step was used, thus requiring 30 time steps to obtain a 1.5-hour fore- 
cast. The method used was the leap-frog scheme described by Haltiner 
(upbl.), with a forward time step used for the first step. To 
eliminate the possibility of obtaining two separate solutions by the 
leap-frog method, the two solutions were averaged every 10 time steps 
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(30 minutes) and then one pass through the averaged field was made 
utilizing a Laplacian smoother with a .04 coefficient for smoothing. 
The leap-frog forecast me thod was then continued. 

The forecasting of temperature and relative humidity by the simple 
advection schemes described in Section 3C does not involve any 
iterative type computations since the wind components were held 
constant for the forecast interval. Consequently, a single 90-minute 
forward time step was used. However, future modifications, which 
should include provisions for obtaining rotational and divergent com- 
ponents of the wind, should also update the wind field with each time 
step so that forecast values of the wind field could be used for 
temperature and relative humidity advection. 
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4. INITIALIZATION OF THE MODEL 



la order to commence computations with the previously discussed 
model, 35 fields of data are needed (7 levels each of height, tempera- 
ture, relative humidity, wind direction and wind speed). Such a 
volume of data would be difficult to obtain by manual means for 
real-time operational use of a model. Therefore, a more realistic 
approach is to use an objective analysis scheme by the computer. 

Although many complicated and very sophisticated analysis schemes 
can be and have been developed (for example see Hughes, 1967), this 
was not the purpose of the research reported in this paper. Con- 
sequently, a simple objective analysis scheme was devised which uses 
the Laplacian operator and an eight-point averaging technique in 
alternating steps. Comparisons of computer analyses and hand analyses 
showed the patterns to be very similar. The major difference between 
the two analysis methods was the gradient in the vicinity of closed 
centers. The hand analyses tended to subjectively spread the gradient 
out over a much larger area whereas this was not so with the computer 
analyses. Which analysis is more nearly correct is merely a matter 
of subjective judgement. With only nine data points within a 14,400 
square mile area during conditions of severe storms, this question 
will probably not be answered until the data network becomes more 
detailed. Nevertheless, the computer analyses were judged to yield 
adequate data fields for testing of the computer program. 

In order to start the analysis a first guess was needed at each 
grid point. This was accomplished by dividing the grid into nine 
"regions of influence" (see Figure 3), each region containing one 
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Figure 3: Regions of Influence. 
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of the nine radiosonde stations. The value of a given parameter at 
each grid point within a region was then set equal to the known value 
at the radiosonde station. 

With the initial guess completed, the analysis was started. The 
first step was to alter the value at the grid points (except for the 
nine data points) by setting the value at that point equal to the 
average of the surrounding eight points according to the following 
diagram : 

4 3 2 

• • • 

5 0 1 

• • • 

6 7 8 

?1 + P 2 + P 3 + P 4 + P 5 + P 6 + P ? + P g (4 -l) 

P 0 8 

After making one pass through the grid in this manner, the Laplacian- 
type scheme was used such that 

P 0 = P„ + 6 (Pi + p 3 + p 5 + P ? - 4P 0 ), (4-2) 

where the value of 6 = 0.1 was arbitrarily selected. After making 
ten passes through the grid in this manner, the eight-point averaging 
technique was again employed. The entire procedure was repeated in 
such a manner that the eight-point average was employed 11 times while 
the Laplacian scheme was utilized 100 times. The resulting analysis, 
after boundary conditions were set, was the analysis used in the pre- 
diction model. The boundary condition imposed was that the value at 
a grid point on an outside row or column was equal to the value at the 
adjacent grid point on the next row or column toward the center of the 
grid. 
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Figures 4a and 5a show hand analyses for two different levels and 
Figures 4b and 5b show the corresponding computer objective analyses. 

A copy of the program for the objective analysis scheme is included in 
Appendix C. 
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Figure 4a: Hand analysis for 950 mb at 1700Z on 30 May 1967. 

Contours every 5 meters. 
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Figure 4b: Computer objective analysis for 95( mb at 1700Z on 

30 May 1967. Contours every 5 meters. 
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Figure 5a: Hand analysis for 500 mb at 1700Z on 30 May 1967. 

Contours every 10 meters. 
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Figure 5b: Computer objective analysis for 500 mb at 1700Z on 

30 May 1967. Contours every 10 meters. 
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5. RESULTS OF THE COMPUTATIONS 



Two sets of data have been utilized with the previously described 
model. The times of the two sets of data are 1700Z, 30 May 1967 and 
2 130Z 5 30 May 1967. 

In synoptic-scale models, the terms in equation (3-la) involving 
vertical motion and divergence have often been neglected since they 
may be considered relatively unimportant by scale analysis (Haltiner, 
upbl.) in the large scale motions of the atmosphere. However, when 
dealing with cumulus -sea le data, this may not be true. In order to 
evaluate the relative importance of the various terms in the model 
throughout the entire grid, computations of the orders of magnitude 
of the various individual terms of equation (3-3) have been performed 
for both the 2130Z and 1700Z sets of observed data for the seven levels 
in the vertical described in Section 1. Tables 2 and 3 show the mini- 
mum and maximum values of each term based on the observed data, not 
on the forecast values. Although these tables only show the minimum 
and maximum values, they do point out the fact that the divergence 
term is very important. There were many instances when all five 
terms were of the same order of magnitude, that is, when each term 
was just as important as each of the other terms. Table 4 shows a 
sample of the values of each of the five terms as extracted from the 
computer printouts . 

^ . r / r,™ -i dv ch) du £ . 

and were of the same sign 

2584 times, or about 54 $ of the time. Therefore, in the cases tested, 

they were additive more often than not which is further justification 

that they should be retained for small-scale data. 



uuc or 40uu values cescea, 



dx dp 
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Table 2: Minimum and maximum values of components of equation (3-3) for 1700Z data for all 

seven isobaric levels. 
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Table 3: Minimum and maximum values of components of equation (3-3) for 2130Z data for all 

seven isobaric levels. 
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Table 4: Sample of values of the components of equation (3-3) at 950 mb for 1700Z data. 



Vertical motion computations yielded considerably different values 
for the 1700Z data and the 2130Z data. Table 5 shows the maximum up- 
ward and downward vertical motions for each isobaric level at each of 
the times. Although the values in Table 5 may, at first, seem somewhat 
large, a transformation by the formula 

co ^ -p gw (5- 1) 

where p is replaced by p/RT from the equation of state to yield 

“ ^ " RT gW ’ < 5 ' 2 ) 

shows that a value of a) = 200 mb hr ^ upward corresponds to a value 

of w ~ .86 m sec ^ upward (at 500 mb) which is certainly a reasonable 

value for cumulus convection. Although values this large were only 

present at a few points, further research is needed to determine what 

reasonable values are for this scale of motion. In general, the 

vertical motion fields showed fairly strong upward vertical motions 

throughout the grid at 1700Z and correlated fairly well with the 

squall- line that was present at that time. There were a few areas 

in which the correlation was not too good but this may possibly be 

due to the fact that the nine radiosondes were not all launched 

simultaneously and a difference of 15 minutes between launch times 

could possibly result in the bad correlation regions. Although the 

radar pictures indicated that dissipation of the squall-line had 

begun at 2130Z by the fact that isolated echoes were present rather 

than a large mass of clouds, the vertical motion fields did not show 

extensive regions of downward motions except at 300 mb. The 2130Z 

values were, however, considerably less in intensity than they were 

at 1700Z which might have been anticipated from the results shown in 

Table 5. 
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Table 5: Maximum upward and downward vertical motions for each isobaric level at 1700Z 

and 2130Z on 30 May 1967 in mb hr"^. 



The significance of the much larger values of vertical motion at 
1700Z as compared to 2130Z was not investigated thoroughly, however 
larger values of divergence were evident (as would be expected) at 
1700Z. This was particularly true at 200 mb, and this would tend to 
be propagated downward to 300 mb and 500 mb by the computation procedure. 

The simple method of forecasting temperature resulted in values 
which verified slightly better than a persistence forecast. The 
average absolute error (not considering whether the error was an 
overestimate or an underestimate) for 97 persistence temperature 
forecasts was 1.8C while for 105 50$ temperature advection forecasts 
the average absolute error was 1.4C. The reason for the different 
number of verifications is due to the fact that some of the 1700Z and 
2130Z radiosonde data terminated before reaching 150 mb. Verifications 
were performed only at the nine radiosonde stations, and only at each 
of the seven discrete isobaric levels. 

The simple method of forecasting relative humidity resulted in 
values which verified essentially the same as a persistence forecast. 

The average absolute error for 97 persistence relative humidity fore- 
casts was 18.0$ while for 105 50$ relative humidity advection forecasts 
the average absolute error was 18.2$. 

The results of the height forecasts did not verify as well as was 
expected. This is undoubtedly due in part to the fact that additional 
refinements are necessary in the prediction method to obtain a con- 
sistent set of equations. Table 6 compares the average absolute 
errors of persistence and model 1.5-hour forecasts for both the 1700Z 
and the 2130Z sets of data. As can be seen in the table, the model 
did perform slightly better than persistence 50$ of the time. 
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AVERAGE ABSOLUTE ERROR (METERS) 
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Table 6: Average absolute errors of height fields (meters) for persistence and model at 

1700Z and 2 130Z for all seven isobaric levels. 



As previously mentioned, the W.D. model yields profiles of vertical 
velocity, hydrometeor water content and temperature excess. In addition, 
total rain expected from the cloud, duration of the precipitation, up- 
draft area, and a profile of updraft radius are included in the results 
of the computations. Figure 6 shows a sample of the graphical output 
of the W.D. model. These results are printed out in graphical form 
by the two subroutines (GRPHCL and GRAPH) which were included with the 
W.D. model (Weinstein and Davis, 1967). The horizontal axis (shown in 
Figure 6 at 1/2 scale) shows the scale to be used for the various para- 
meters and the vertical axis on the left shows the corresponding 
heights in meters, with the bottom figure (in this case 1642) being 
the first increment (200 meters) above cloud base. In the actual 
computer output of Figure 6, the lines are represented by letters, 
such that the vertical velocity curve is a series of the letter ,f W ff , 
hydrometeor water content of the letter M Q M , temperature excess of 
the letter "T", updraft radius of the letter "R", and the zero axis 
for the temperature excess curve is represented by the symbol ff $ M . 

The W.D. model has a provision for setting a constant value for 
the updraft radius or for allowing it to vary within a cloud. How- 
ever, for this research the updraft radius was always considered to 
be constant and equal to 1 km. The profile shown in Figure 6 was 
computed based upon the 1830Z, 30 May 1967 forecasted values of 
temperature, relative humidity and height of the cloud base, and 
by assuming that cloud base was at 850 mb for station CHK. 

Many soundings were used as input to the W.D. model, and in many 
cases the values of the temperature excess determined by the model 
seemed to be too large. Values as high as 4 or 5 degrees Celsius 
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occurred quite often, with an occasional value of 7 or 8 degrees 
Celsius. This intuitively seems to be too high but additional 
observations are needed to determine if such values actually exist 
within cumulus clouds. In addition, apparent discontinuities in 
temperature excess, as evidenced near 8000 meters on Figure 6, are 
quite common in the graphical output of the model. These large tem- 
perature excess values and temperature excess discontinuities occurred 
not only with forecast soundings but also with actual observed soundings 
from the NSSL data. 

Figure 6 is based upon an ice nucleation temperature of -25C which 
is assumed to be representative of natural conditions in the atmos- 
phere. The W.D. model allows for changing the ice nucleation 
temperature to a value of, say, -6C to represent the conditions that 
might occur if cloud seeding were to be undertaken. However, this 
is presently not included in the subroutine version of the W.D. 
model (subroutine CLOUDS) since the reason for using it as part of 
a prediction model was to aid in severe storm forecasting. 
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6. SUGGESTIONS FOR FUTURE DEVELOPMENTS AND IMPROVEMENTS 



The first modification to the previously discussed model, which 
became evident during testing of the model, would be to expand the 
horizontal grid dimensions to at least 32 x 32 . Figure 3 illustrates 
that four of the nine grid points lie within the first two rows of 
the boundaries. This more or less negates application of the prognostic 
equations to these grid points due to finite differencing. It is esti- 
mated that this improvement would increase the storage requirements of 
the program by about 50$, and since the program is already quite large, 
modifications would have to be made to allow for the sharing of some 
storage locations wherever possible. 

At present the model lacks the proper provisions for predicting 
the complete wind field and for updating it after each time step. 

A consistent scheme for updating the wind field should be included 
such that the predicted winds could be recovered every time step. 

In addition, this would allow predictions further ahead in time as 
well as allow for updating the vertical motion fields with new 
divergence fields . 

As discussed in Section 3A, the present version of the forecast 
model does not provide for vertical coupling of the seven isobaric 
levels since the vorticity in the vertical advection term and the u 
and v components of the wind in the twisting term were held constant 
for the 1.5-hour forecast interval. To correct this deficiency would 
involve performing one time step at a given level before proceding to 
the next level. Also, additional storage would be required to save the 
field which operates on for one time step back so that vertical 
derivatives could be taken at each level on fields which all correspond 
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to the same time. As far as programming considerations are concerned, 
this essentially involves nesting the vertical level DO LOOP inside 
of the time step DO LOOP. 

When the model has been modified to allow for vertical coupling of 
the seven isobaric levels, additional levels in the vertical could be 
included so that vertical layer depths would not exceed 100 mb. This 
would require perhaps 10 or 11 levels which would increase the storage 
requirements of the computer considerably, however, this would allow 
a much finer detail for the forecasted soundings. In addition, it 
should aid in the delineation of the vertical velocities at each 
leve 1 . 

A consistent dynamical method should be included for the prediction 
of temperature and relative humidity. The present method is more 
empirical than dynamical and better verifications should be obtained 
if a more realistic approach is devised. As an alternative, a better 
constant (vice 50 fo) might be tested, particularly in the case of 
relative humidity where the errors were negative much more often 
than they were positive. 

In subroutine SETUP (see Appendix C), cloud base is assumed to be 
at 850 mb for entry into subroutine CLOUDS (see Appendix C). A routine 
could be developed to be included in SETUP such that the lifting con- 
densation level, mixing condensation level or convective condensation 
level would be computed as the cloud base before entering the W.D. 
model (subroutine CLOUDS). Comparison of results obtained for cloud 
base at 950 mb and 850 mb indicate that the cloud base and the associ- 
ated temperature and relative humidity are very critical to the results 
obtained from the W.D. model. 
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Once a complete internally consistent model is set up, calculations 
out past 1.5 hours could be tried to determine at what point such a 
model loses its skill. A few runs were made out to 4.5 hours and the 
results appeared to indicate values of height which were too large, and 
the temperature fields appeared to become unstable. 

Kessler (1967) proposed the formula for the fall speed of pre- 
cipitation to be of the form 

V = 5.0939 M " 125 . (6-1) 



The corresponding formula in the W.D. model is 

125 

V = 15.39 Q u . 

n 



(6-2) 



The units of M are gm m whereas the units of are gm gm . In 
order to convert M to the units of Q , it is necessary to divide by 
the density of air, p. Weinstein and Davis (1967) have apparently 
selected a representative value for p. It is recommended that future 
modifications of the W.D. model include the density in the computation. 
Since M = pQ^ , and p - p/RT from the equation of state, it would be a 
simple change in the programming to use the following formula for the 
hydrometeor terminal velocity: 



V = 5.0939 




(6-3) 



Since the duration calculation is based on the hydrometeor terminal 
velocity at cloud top, (6-3) will permit values of V which not only 
depend on the hydrometeor liquid water content but also will depend 
on the pressure and temperature at the cloud top. 

Finally, it has been suggested to the author by Dr. J. D. Mahlman 
(personal communication) that the only solution to overcome the many 
difficulties with this scale of motion would be a direct integration 
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of the primitive equations. It is believed that no such research has 
been attempted for mesoscale data and this approach definitely warrants 
consideration for future attempts at numerical modelling of such scales 
of motion. 
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7 . SUMMARY 



A numerical weather prediction model has been described for use 
with cumulus-scale data. Together with the W.D. model, it provides 
forecasts of height, temperature, and relative humidity at seven 
isobaric levels as well as vertical profiles of vertical velocity, 
temperature excess and hydrometeor liquid water content at 1.5-hour 
intervals. The model was tested on two sets of data and yielded fore- 
casts which were, in general, statistically about the same as a per- 
sistence forecast. This can probably be attributed in part to the 
fact that the present version has no provisions for coupling the 
seven isobaric levels nor for updating the vertical motion fields 
after each time step. Another factor is that the temperature and 
relative humidity forecast schemes were more empirical than dynamical, 
and were not consistent witn the vorticity equation. This may account 
for the lack of skill in forecasting these parameters. Suggestions 
were made to correct some of these inconsistencies in the event that 
further research is conducted with the model. However, the possibility 
of applying numerical forecasting techniques to small-scale data does 
appear to have possibi ities but extensive research is needed in this 
area before useful forecasts will be obtained. 

The model requires approximately 5.5 minutes on the IBM OS/360 
(MVT versions i5',- 16) computer to yield a 1.5 hour forecast. An 
additional 14 seconds per field is required for the objective analysis 
scheme described in Section 4, thus requiring about 13.5 minutes for 
the analysis of 35 fields of data and the prediction of the afore- 
mentioned parameters for 1.5 hours. 
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APPENDIX A 



THE NSSL RADIOSONDE NETWORK 

There were 10 radiosonde stations within the NSSL region considered 
for this study which are listed below. The grid-point numbers associ- 
ated with each station are shown on Figure 3, where the stations are 
shown by their call letters. 



Call letters 


Name 


Location 
Lat. Long. 


SPS 


Sheppard AFB , Texas 


33.97 


98.48 


RIN 


Ringling, Oklahoma 


34.17 


97.58 


LTS 


Altus, Oklahoma 


34.70 


99.33 


FSI 


Fort Sill, Oklahoma 


34.65 


98.40 


PVY 


Pauls Valley, Oklahoma 


34.70 


97.22 


COR 


Cordell, Oklahoma 


35.30 


98.97 


CHK 


Chickasha , Ok lahoma 


35.10 


97.97 


OKC 


Oklahoma City, Oklahoma 


35.40 


97.60 


TIK 


Tinker AFB, Oklahoma 


35.42 


97.38 


WAT 


Watonga, Oklahoma 


35.85 


98.42 



All of these stations, except OKC , were used in this study. The 
reason for not including data from OKC is that their data corresponds 
to synoptic data times rather than at 1.5-hour intervals. 
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APPENDIX B 



ON THE ELLIPTICITY REQUIREMENT OF THE BALANCE EQUATION 



The ellipticity requirement of the balance equation (Charney, 
and Bolin, 1955) was given by equation (3-5) which is 




1955 



(3-5) 



Bolin (1955) stated that the data at grid points which did not satisfy 
(3-5) was modified so that (3-5) was satisfied. This section investi- 
gates the extent to which grid distance influences the requirement of 
(3-5). 

Expanding (3-5) results in 

^ - V 2 z + y 2 > 0 (B-l) 

d f 

m 

where 

2 

V = finite difference form of the Laplacian operator, 
m = map factor^ 1.19, 
g = gravity « 9.80, 
d = grid distance, 
z = height , 

-4 

and f = mean Coriolis parameter = .836x10 
m 

Upon substitution of the above values for each factor, (B-l) reduces to 

V 2 z + .00087d 2 > 0 (B-2) 

where z is in meters and d is in nautical miles. Substituting a 
value of d = 5 nm yields 

V 2 z > -.02175. (B-3) 

By measuring heights to the nearest meter, the only way to satisfy 

2 

(B-3) is to have V z > 0 . To satisfy (B-3) on a grid with d = 5 nm 
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requires z to be measured at least to the nearest .01 m, and this would 

still not guarantee meeting the requirement. It is easily shown with 

equation (B-2) that the smallest grid distance on which (3-5) can be 

2 

satisfied, without eliminating the possibility of V z < 0 and still 
measure z to the nearest meter, is 35 nm. 

The result from above is what prompted the substitution of the 
measured wind field into the balance equation in order to force it 
to be elliptic. Tests were made to determine with what frequency the 
balance equation ellipticity requirement was satisfied with this scale 
of data. The results show that out of 6776 points, only 62.4 $ satis- 
fied the requirement. At 1700Z the requirement was satisfied 59.4$ 
of the time while at 2130Z the figure was 65.5$. 
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APPENDIX C 






This appendix includes a listing of the computer program used for 
the research described by the preceding sections . In order to aid 
any possible research in the future with this program, each subroutine 
will be listed with a brief explanation of the calling arguments. 

1. Subroutine VERTMO 

U = Zonal component of the horizontal wind in m sec ^ . 

V = Meridional component of the horizontal wind in m sec 
OMEGA = O) = Vertical component of the wind in (x, y, p, t) 
coordinate system in mb sec 
NI = Number of grid points along x-axis. 

NJ = Number of grid points along y-axis. 

NK = Number of levels in the vertical. 

DIST = Distance between grid points in meters. 

DSI = The negative of the divergence field for 1000 mb = the 
sum of the divergence from 950 mb to 150 mb. 

2. Subroutine FORFUN 

VORT = Vertical component of relative vorticity. 

F = Field for forcing function. 

FMAP = Map factor (mean). 

FBAR = Coriolis parameter (mean). 

DIST = Distance between grid points in meters. 

DT = Forward time step in seconds. 

USPD = Zonal component of the horizontal wind in m sec 
VSPD = Meridional component of the horizontal wind in m sec 
K = The level the main part of the program is at when this 
subroutine is called. 
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IZ = Counter to determine at what time step (1-30) the main 
program is . 

2 -1 

PSI = The i|f field in m sec 

OMEGA = cd = Vertical component of the wind in (x, y, p, t) 
coordinate system in mb sec 
3 . Subroutine RELAXI 

M = Number from 2, NI = a number one greater than what is 
needed in DO LOOP. 

N = Same as M except 2, NJ. 

F = At field for relaxation of balance equation. 

A = t field in m^ sec 
B = Z field in meters. 

L = Counter to keep track of number of points which have 
converged . 

EPS = Epsilon = convergence criteria. 

ALPHA = Over-relaxation constant. 

FMAP = Map factor (mean). 

FBAR - CoriolLs parameter (mean). 

DIST = Distance between grid points in meters. 

IL ™ Lower point for DO LOOP in x direction. 

JL = Lower point for DO LOOP in y direction. 

LM = Iteration number from main program. 

U = Zonal component of the horizontal wind in m sec 
V = Meridional component of the horizontal wind in m sec ^ . 
K = Level in the vertical from main program. 

G = Gravity. 
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4. Subroutine RELAX 



Same as in RELAXI except : 

F = Forcing function. 

R = Residual field. 

ICOUNT = L. 

II = IL. 

Jl = JL. 

5. Subroutine SMOOTH 

D = Field to be smoothed. 

\ 

RF = Smoothing factor (<1.0). 

IPASS = Number of smoothing passes to make through the grid. 
IL, JL = Lower point for DO LOOP. 

IH, JH = Upper point for DO LOOP. 

6. Subroutine PROG 
TEMP = Temperature. 

RH = Relative humidity. 

DIR = Wind directions. 

SPD = Wind speeds . 

USPD = Zonal component of SPD in m sec 

VSPD = Meridional component of SPD in m sec ^ . 

DIST = Distance between grid points in meters. 

NK = Number of levels in the vertical. 

NI = Number of grid points on x-axis . 

NJ = Number of grid points on y-axis . 

Nil = NI-1. 

NJ1 = NJ-1. 

OMEGA = oo. 

Z = Height fie Id . 
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7. Subroutine SETUP 



TEMP = Temperature. 

RH = Relative humidity. 

Z = Height fie Id . 

IL = A number from 1 to NI, usually 1. 

IH = A number from 1 to NI , usually NI . 

JL = A number from 1 to NJ, usually 1. 

JH = A number from 1 to NJ, usually NJ. 

8. Subroutine CLOUDS. 

This is the modified version of the W.D. model (Weinstein and 
Davis, 1967). 

PDATA = Pressure levels. 

ZDATA = Height field. 

RHDATA = Relative humidity field. 

ZZ1 = Height of cloud base. 

DZ = Vertical increment (meters) to be used. 

NCALL = Number of times this subroutine called. Used as 
variable to bypass a READ statement. 

IPT = I coordinate for set of data coming in. 

JPT = J coordinate for set of data coming in. 

9. Calling arguments of GRPHCL and GRAPH are all parameters which 
have been computed or set in subroutine CLOUDS, 

10. Subroutine METMAP 

This subroutine was not written by the author, but is part 
of the Naval Postgraduate School Computer Facility Library and has 
been included for information. It is simply a shading routine that 

will print out a field of data with a 0.5 inch grid spacing and contour 

the field. 
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Y = Two-dimensional field to be contoured. 

N = Number of rows I in the array to be contoured. 

M = Number of columns J in the array to be contoured. 

T = Title for printout . Up to 96 columns of alpha-numeric 
information . 

BND = Bandwidth desired for contouring. 

AZ = Scaling constant; that is, each value will be multiplied 

by AZ before printing it out. Only the first three numbers 
to the right of the decimal point are on the output, so AZ 
can be used to control which numbers are to be printed out. 

BZ = Additive constant. This is useful when working with D-values. 
Usually, however, it is 0.0. If used, this constant will be 
added to each value of Y before printing it out. 

AMIN = Minimum value for subroutine to start contouring. 

IJT = 0 means subroutine will compute the minimum value and 
start contouring at that value. A value must still be 
specified for AMIN but it will not be used. 

ICON = 1 if contouring desired. 

= 0 if no contouring desired. 

When working with this subroutine, it must be remembered that the 
grid point I = 1, J = 1, is the upper left corner of the grid, and I 
increases downward while J increases to the right. If an array is 
defined with the lower left corner as the point (1, 1), the follow- 
ing set of FORTRAN statements will get the array in the proper order 
for the subroutine: 
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DO 5 1=1, 24 
DO 5 J=l, 24 

5 DUMMY (J, I)=Z (I ,24- J+l) 

The field called DUMMY is then taken into METMAP as Y. 
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